#include "get_arg.h"
#include <boost/filesystem/path.hpp>
#include <boost/filesystem/operations.hpp>
#include <iostream>
#include <stdlib.h>
#include <sstream>
#include <fstream>
#include <matrix.h>
#include <matrix_io.h>
#include <portable_timer.h>
#include <boost/lexical_cast.hpp>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <tree_UPGMA.h>
#include <tree.h>


using namespace std;
using namespace bcpp_matrix;
using boost::lexical_cast;
using std::ostringstream;
using std::istringstream;
using std::ofstream;
using std::ifstream;

unsigned get_file_size(const string& name)
{
	struct stat file_stat;
	if(stat(name.c_str(), &file_stat))
		return 0;
	return file_stat.st_size;

}

struct t_iter_index
{
	unsigned index;
	vector<string>::const_iterator iter;
	t_iter_index(unsigned i, const vector<string>::const_iterator& it)
		:index(i), iter(it){}
	bool operator<( const t_iter_index& other)
	{
		return *(other.iter) < *iter;
	}
	operator unsigned () const {return index;}

};

void sort_matrix_by_name(t_matrix_f& matrix)
{
	vector<t_iter_index> named_indices;
	for (unsigned i = 0; i < matrix.names.size(); ++i)
		named_indices.push_back(t_iter_index(i, matrix.names.begin() + i));
	sort(named_indices.begin(),  named_indices.end());

	vector<unsigned> indices(named_indices.begin(),  named_indices.end());
	matrix.reduce_to_subset(indices);

}

struct t_output_message_with_timer
{
	t_portable_timer timer;
	string end_message;
	unsigned size;
	unsigned pos;
	void set_size(unsigned i) {size = i;}
	t_output_message_with_timer(const string& message, bool delay = false)
		: size(0), pos(0)
	{
		if (delay)
		{
			end_message = message;
		}
		else
		{
			pos = message.length() + 8;
			cerr << "\t" << message;
		}
	}
	~t_output_message_with_timer()
	{
		if (end_message.length())
		{
			pos += end_message.length() + 8;
			cerr << "\t" << end_message;
		}

		if (size)
		{
			ostringstream str;
			if (size > 1024*1024*10)
				str << " (" << size / 1024 / 1024 << " MB) ";
			else if (size > 1024 * 10)
				str << " (" << size / 1024 << " kB) ";
			else
				str << "(" << size << " bytes) ";
			pos += str.str().length();
			cerr << str.str();
		}
		if (pos < 70)
			cerr << string(70 - pos, ' ');
		cerr << timer.pretty_reset() << "\n";
	}
};


void output_matrix_bool(const t_dense_matrix_bool& m, ostream& os)
{
	for (unsigned i = 0; i < m.size(); ++i)
	{
		for (unsigned j = 0; j < m.size(); ++j)
			os << (m[i][j] ? "1 " : "0 ");
		os << '\n';
	}
}

void test_tree_read(std::istream& f_temp, const std::string& tree_str, unsigned& cnt_lines)
{
	bcpp_tree::t_tree input_tree;
	input_tree.set_sort_by_names(true);
	std::string input_tree_str;
	std::string tree_name;
	bool success = read_tree_from_stream(f_temp, input_tree, cnt_lines);
	assert(success);
	input_tree_str = ">" + input_tree.get_name() + "\n" + input_tree.to_nh_str();
	if (input_tree_str != tree_str)
	{
		cerr << "\n>>original tree_str = [" << tree_str << "]\n";
		cerr << "\n>>read tree_str = [" << input_tree_str << "]\n";
	}
	assert(input_tree_str == tree_str);

}
int main (int argc, char *argv[])
{
	t_program_args args;
	if (!process_command_line_options(argc, argv, args))
		return 1;

	const unsigned MATRIX_SZ = args.size;
	const float DEFAULT_VAL = args.default_val;
	const unsigned cnt_values = MATRIX_SZ * (MATRIX_SZ - 1) / 2;
	const unsigned PROPORTION = min(cnt_values, unsigned(cnt_values * args.fill));


	//
	// Allocate elements
	//
	t_matrix_f matrix1;
	{
		cerr << "Allocate " << MATRIX_SZ << " elements\n";
		matrix1.allocate(MATRIX_SZ, DEFAULT_VAL);
		// DEBUGG
		cerr << matrix1.null_value << "--??\n";
	}

	//
	// Assign names
	//
	{
		t_output_message_with_timer message("Assign names");
		for (unsigned i = 0; i < MATRIX_SZ; ++i)
			matrix1.names[i] = "A_" + lexical_cast<string>(i);
		assert(matrix1.matrix.size() == MATRIX_SZ);
	}
	matrix1.matrix_name = "test1";


	//
	// Randomly assign values
	//
	{
		t_output_message_with_timer message("Randomly assign values");

		// diagonals are zeros to be compatible with triangle format
		for (unsigned ii = 0; ii < MATRIX_SZ; ++ii)
		{
			matrix1.matrix[ii][ii] = 0.0;
		}

		// first make sure no columns are entirely null
		std::vector<bool> used_names(MATRIX_SZ, false);
		unsigned cnt_filled = MATRIX_SZ;
		for (unsigned ii = 0; ii < MATRIX_SZ; ++ii)
		{
			while (1)
			{
				// get random other column
				unsigned jj = rand() %  MATRIX_SZ;
				if (ii == jj )
					continue;
				matrix1.matrix[ii][jj] =
				matrix1.matrix[jj][ii] = rand() * 10.0 /  RAND_MAX;
				++cnt_filled;
				break;
			}
		}

		// keep going until PROPORTION cells filled
		while (cnt_filled < PROPORTION)
		{
			unsigned ii = rand() %  MATRIX_SZ;
			unsigned jj = rand() %  MATRIX_SZ;
			if (ii == jj)
				continue;
			if (!matrix1.is_missing(ii, jj ))
				continue;
			matrix1.matrix[ii][jj] =
			matrix1.matrix[jj][ii] = rand() * 10.0 /  RAND_MAX;
			++cnt_filled;
		}


	}


	matrix1.count_missing();
	cerr << "\t\t" << matrix1.get_count_missing() << " missing out of ";
	cerr << matrix1.get_count_pairs() << " expected\n";
	// DEBUGG
	cerr << matrix1.null_value << "--??\n";

	//
	// Sort names
	//
	{
		t_output_message_with_timer message("Sort names");
		sort_matrix_by_name(matrix1);
	}

	//
	// Round values so that comparisons will work
	//
	{
		t_output_message_with_timer message("Round up values by writing out then reading back");
		stringstream strm;
		write_matrix(strm, matrix1, eMIO_SQUARE, DEFAULT_VAL);
		read_matrix(strm, matrix1, eMIO_SQUARE, DEFAULT_VAL);
		{
			ofstream f_temp("DEBUG");
			write_matrix(f_temp, matrix1, eMIO_SQUARE, DEFAULT_VAL);
		}
	}

//	write_matrix_to_square(cout, matrix1, true, false);

	// DEBUGG
	cerr << matrix1.null_value << "--??\n";
	//
	// binary format
	//
	{
		cerr << "\n";
		t_output_message_with_timer message("Matrix in binary format", true);
		string data;
		t_matrix_f matrix2;
		{
			ofstream f_temp(args.temp_file.c_str());
			if (!f_temp)
			{
				cerr << "Error:\n\tCould not open " << args.temp_file << "\n";
				return 1;
			}
			t_output_message_with_timer message("Write matrix in binary format");
			write_matrix(f_temp, matrix1, eMIO_BINARY, DEFAULT_VAL);
			// DEBUGG
			cerr << matrix1.null_value << "--??\n";
			//write_matrix(f_temp, matrix1, eMIO_SQUARE, false);
			f_temp.close();
			message.set_size(get_file_size(args.temp_file));
		}
		{
			ifstream f_temp(args.temp_file.c_str());
			if (!f_temp)
			{
				cerr << "Error:\n\tCould not open " << args.temp_file << "\n";
				return 1;
			}
			t_output_message_with_timer message("Read matrix in binary format");
			read_matrix(f_temp, matrix2, eMIO_BINARY, DEFAULT_VAL);
		}
		assert(matrix2.names == matrix1.names);
		//for (unsigned i = 0; i < MATRIX_SZ; ++i)
		//{
		//	for (unsigned j = 0; j < MATRIX_SZ; ++j)
		//	{
		//		cerr << "[" << i << "," << j << "] = " << matrix1.matrix[i][j] << "\t"
		//				 << matrix2.matrix[i][j]
		//				<<	(matrix1.is_missing(i, j) ? " missing1": " present1") << "\t"
		//				<<	(matrix2.is_missing(i, j) ? " missing2": " present2") << "\n";
		//	}
		//}
		assert(matrix2.matrix == matrix1.matrix);
		cerr << "matrix1 missing count = " << matrix1.get_count_missing() << "\n";
		cerr << "matrix1 missing count = " << matrix2.get_count_missing() << "\n";
		assert(matrix2 == matrix1);
	}

	//
	// Lower triangular format
	//
	{
		cerr << "\n";
		t_output_message_with_timer message("Matrix in lower triangular format", true);
		string data;
		t_matrix_f matrix2;
		{
			ofstream f_temp(args.temp_file.c_str());
			if (!f_temp)
			{
				cerr << "Error:\n\tCould not open " << args.temp_file << "\n";
				return 1;
			}
			t_output_message_with_timer message("Write matrix in lower triangular format");
			write_matrix(f_temp, matrix1, eMIO_LOWER, DEFAULT_VAL);
			f_temp.close();
			message.set_size(get_file_size(args.temp_file));
		}
		{
			ifstream f_temp(args.temp_file.c_str());
			if (!f_temp)
			{
				cerr << "Error:\n\tCould not open " << args.temp_file << "\n";
				return 1;
			}
			t_output_message_with_timer message("Read matrix in lower triangular format");
			read_matrix(f_temp, matrix2, eMIO_LOWER, DEFAULT_VAL);
		}
		{
			ofstream f_temp(args.temp_file.c_str(), std::ios::app | std::ios::ate);
			write_matrix(f_temp, matrix1, eMIO_SQUARE, DEFAULT_VAL);
			write_matrix(f_temp, matrix2, eMIO_SQUARE, DEFAULT_VAL);
		}
		assert(matrix2 == matrix1);
	}

	//
	// pair-wise format
	//
	{
		cerr << "\n";
		t_output_message_with_timer message("Matrix in pairwise format", true);
		string data;
		t_matrix_f matrix2;
		{
			ofstream f_temp(args.temp_file.c_str());
			if (!f_temp)
			{
				cerr << "Error:\n\tCould not open " << args.temp_file << "\n";
				return 1;
			}
			t_output_message_with_timer message("Write matrix in pair-wise format");
			write_matrix(f_temp, matrix1, eMIO_PAIRS, DEFAULT_VAL);
			f_temp.close();
			message.set_size(get_file_size(args.temp_file));
		}
		{
			ifstream f_temp(args.temp_file.c_str());
			if (!f_temp)
			{
				cerr << "Error:\n\tCould not open " << args.temp_file << "\n";
				return 1;
			}
			t_output_message_with_timer message("Read matrix in pair-wise format");
			read_matrix(f_temp, matrix2, eMIO_PAIRS, DEFAULT_VAL);
		}
		{
			t_output_message_with_timer message("Sort matrix for comparison");
			sort_matrix_by_name(matrix2);
		}
		assert(matrix2 == matrix1);
	}



	//
	// Square format
	//
	{
		cerr << "\n";
		t_output_message_with_timer message("Matrix in square format", true);
		string data;
		t_matrix_f matrix2;
		{
			ofstream f_temp(args.temp_file.c_str());
			if (!f_temp)
			{
				cerr << "Error:\n\tCould not open " << args.temp_file << "\n";
				return 1;
			}
			t_output_message_with_timer message("Write matrix in square format");
			write_matrix(f_temp, matrix1, eMIO_SQUARE, DEFAULT_VAL);
			f_temp.close();
			message.set_size(get_file_size(args.temp_file));
		}
		{
			ifstream f_temp(args.temp_file.c_str());
			if (!f_temp)
			{
				cerr << "Error:\n\tCould not open " << args.temp_file << "\n";
				return 1;
			}
			t_output_message_with_timer message("Read matrix in square format");
			read_matrix(f_temp, matrix2, eMIO_SQUARE, DEFAULT_VAL);
		}
		assert(matrix2 == matrix1);
	}

	boost::filesystem::remove(boost::filesystem::path(args.temp_file));;

	cerr << "\n";

	cerr << "test operators +-*/ ...\n";
	for (unsigned i = 0; i < 50; ++i)
	{
		cerr << ".";
		const unsigned size = rand() %  MATRIX_SZ;
		t_dense_matrix_f matrix2;
		matrix2.allocate(size + 21);
		matrix2.fill(4.0);

		t_dense_matrix_f matrix1;
		matrix1.allocate(size);
		for (unsigned i = 0; i < size; ++i)
			for (unsigned j = 0; j < size; ++j)
				matrix1[i][j] = rand() % 1000;

		t_dense_matrix_f matrix3 = (matrix1 * matrix2) / matrix2;
		assert(matrix3 == matrix1);
		matrix3 = (matrix1 + matrix2) - matrix2;
		assert(matrix3 == matrix1);
		matrix3 = matrix1;
		matrix3 *= matrix2;
		matrix3 /= matrix2;
		matrix3 += matrix2;
		matrix3 -= matrix2;
		assert(matrix3 == matrix1);
	}
	cerr << "\n";

	cerr << "test operators | & ...\n";
	{
		cerr << ".";
		const unsigned size = 5;
		t_dense_matrix_bool matrix2;
		matrix2.allocate(size);
		for (unsigned i = 0; i < matrix2.size(); ++i)
			for (unsigned j = 0; j < matrix2.size(); ++j)
				matrix2[i][j] = bool (rand() % 2);

		t_dense_matrix_bool matrix1;
		matrix1.allocate(size + 2);
		for (unsigned i = 0; i < matrix1.size(); ++i)
			for (unsigned j = 0; j < matrix1.size(); ++j)
				matrix1[i][j] = bool (rand() % 2);

		t_dense_matrix_bool matrix3 = matrix1 | matrix2;

		cerr << ">matrix1\n";
		output_matrix_bool(matrix1, cerr);
		cerr << ">matrix2\n";
		output_matrix_bool(matrix2, cerr);
		cerr << ">1 | 2\n";
		output_matrix_bool(matrix3, cerr);
		cerr << ">2 | 1\n";
		output_matrix_bool(matrix2 | matrix1, cerr);
		cerr << ">1 & 2\n";
		output_matrix_bool(matrix1 & matrix2, cerr);
		cerr << ">2 & 1\n";
		output_matrix_bool(matrix2 & matrix1, cerr);
	}
	cerr << "\n";





	//
	// Multiple formats
	//
	{
		cerr << "\n";
		t_output_message_with_timer message("Matrix in multiple formats", true);
		string tree_str;
		{
			args.VERBOSE() << "\tBuilding UPGMA tree...\n";
			bcpp_tree::t_tree upgma_tree;
			upgma_tree.set_sort_by_names(true);
			t_progress_indicator dots(args.VERBOSE(), 1);
			upgma_tree::upgma_make_tree(upgma_tree, matrix1.names, matrix1.matrix,
										false /*args.fill_lacunae*/, dots);
			tree_str = ">" + matrix1.matrix_name + "\n" + upgma_tree.to_nh_str();

		}

		string data;
		t_matrix_f matrix2;
		{
			ofstream f_temp(args.temp_file.c_str());
			if (!f_temp)
			{
				cerr << "Error:\n\tCould not open " << args.temp_file << "\n";
				return 1;
			}
			{
				matrix1.sort_by_name_alpha();
				t_output_message_with_timer message("Write matrix in multiple formats");
				cerr << "!1\tSquare format\n";
				write_matrix(f_temp, matrix1, eMIO_SQUARE, DEFAULT_VAL);
				f_temp << tree_str<<"\n";
				cerr << "!1\tBinary format\n";
				write_matrix(f_temp, matrix1, eMIO_BINARY, DEFAULT_VAL);
				f_temp << tree_str<<"\n";
				cerr << "!1\tLower format\n";
				write_matrix(f_temp, matrix1, eMIO_LOWER, DEFAULT_VAL);
				f_temp << tree_str<<"\n";
				cerr << "!1\tPairs format\n";
				write_matrix(f_temp, matrix1, eMIO_PAIRS, DEFAULT_VAL);
				f_temp << tree_str<<"\n";
				cerr << "!2\tSquare format\n";
				write_matrix(f_temp, matrix1, eMIO_SQUARE, DEFAULT_VAL);
				f_temp << tree_str<<"\n";
				cerr << "!2\tBinary format\n";
				write_matrix(f_temp, matrix1, eMIO_BINARY, DEFAULT_VAL);
				f_temp << tree_str;
				cerr << "!2\tLower format\n";
				write_matrix(f_temp, matrix1, eMIO_LOWER, DEFAULT_VAL);
				f_temp << tree_str;
				cerr << "!2\tPairs format\n";
				write_matrix(f_temp, matrix1, eMIO_PAIRS, DEFAULT_VAL);
				f_temp << tree_str;
				cerr << "!3\tSquare format\n";
				write_matrix(f_temp, matrix1, eMIO_SQUARE, DEFAULT_VAL);
				f_temp << tree_str;
			}
			f_temp.close();
			message.set_size(get_file_size(args.temp_file));
		}
		{
			ifstream f_temp(args.temp_file.c_str());
			if (!f_temp)
			{
				cerr << "Error:\n\tCould not open " << args.temp_file << "\n";
				return 1;
			}
			{
				t_output_message_with_timer message("Read matrix in multiple formats");
				std::ostringstream ignore;
				t_progress_indicator dots(ignore, 1);
				unsigned cnt_lines = 0;

				cerr << ">1\tSquare format\n";
				f_temp >> t_matrix_istrm_wrapper_f(matrix2, eMIO_SQUARE,  cnt_lines,  dots);
				//read_matrix(f_temp, matrix2, eMIO_SQUARE, ignore_cnt_duplicates, DEFAULT_VAL);
				assert(matrix2 == matrix1);
				test_tree_read(f_temp, tree_str, cnt_lines);

				cerr << ">1\tBinary format\n";
				read_matrix(f_temp, matrix2, eMIO_BINARY, DEFAULT_VAL);
				assert(matrix2 == matrix1);
				test_tree_read(f_temp, tree_str, cnt_lines);

				cerr << ">1\tLower format\n";
				read_matrix(f_temp, matrix2, eMIO_LOWER, DEFAULT_VAL);
				assert(matrix2 == matrix1);
				test_tree_read(f_temp, tree_str, cnt_lines);

				cerr << ">1\tPairs format\n";
				read_matrix(f_temp, matrix2, eMIO_PAIRS, DEFAULT_VAL);
				matrix2.sort_by_name_alpha();
				if (!(matrix2 == matrix1))
				{
					std::cerr << "first\n";
					write_matrix(std::cerr, matrix1, eMIO_PAIRS, DEFAULT_VAL);
					std::cerr << "second\n";
					write_matrix(std::cerr, matrix2, eMIO_PAIRS, DEFAULT_VAL);
				}
				assert(matrix2 == matrix1);
				test_tree_read(f_temp, tree_str, cnt_lines);

				cerr << ">2\tSquare format\n";
				read_matrix(f_temp, matrix2, eMIO_SQUARE, DEFAULT_VAL);
				assert(matrix2 == matrix1);
				test_tree_read(f_temp, tree_str, cnt_lines);

				cerr << ">2\tBinary format\n";
				read_matrix(f_temp, matrix2, eMIO_BINARY, DEFAULT_VAL);
				assert(matrix2 == matrix1);
				test_tree_read(f_temp, tree_str, cnt_lines);

				cerr << ">2\tLower format\n";
				read_matrix(f_temp, matrix2, eMIO_LOWER, DEFAULT_VAL);
				assert(matrix2 == matrix1);
				test_tree_read(f_temp, tree_str, cnt_lines);

				cerr << ">2\tPairs format\n";
				read_matrix(f_temp, matrix2, eMIO_PAIRS, DEFAULT_VAL);
				matrix2.sort_by_name_alpha();
				assert(matrix2 == matrix1);
				test_tree_read(f_temp, tree_str, cnt_lines);

				cerr << ">3\tSquare format\n";
				read_matrix(f_temp, matrix2, eMIO_SQUARE, DEFAULT_VAL);
				assert(matrix2 == matrix1);
				test_tree_read(f_temp, tree_str, cnt_lines);
			}
		}
	}




	cerr << "\n";
	cerr << string(10, '!') << "\nEverything worked!!\n";

	return 0;
}


